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1  INTRODUCTION 


Background 

Because  geographic  information  systems  (GISs)  have  been  enhanced  with  more  robust  query 
capabilities,  they  have  evolved  from  basic  spatial  data  storage,  retrieval,  and  display  systems  to 
more  powerful  spatial  modeling  systems.  Languages  that  empower  an  end-user  to  design  and  create 
spatial  models  will  advance  this  development  even  further.  Tomlin  (1990)  described  a  spatial  mod¬ 
eling  language  for  raster  data  with  three  major  classes  of  map  operations:  functions  of  individual 
cells  (e.g.,  map  recoding  and  Boolean  overlay);  functions  of  cells  within  neighborhoods  (e.g.,  filters 
and  slope-aspect  generation);  and  functions  of  cells  within  regions  or  zones  (e.g.,  area  calculations). 
This  classification  could  be  applied  to  image-processing  operations  as  well  because  they  are  inher¬ 
ently  raster- based  and  employ  analysis  methods  very  similar  to  some  GIS  operations  (Burrough 
1986). 

A  GIS  map  algebra  was  needed  that  would  accommodate  Tomlin’s  map  classifications,  al¬ 
low  image-processing  operations  and  could  be  implemented  in  the  Geographical  Resources  Analysis 
Support  System  (GRASS).  GRASS  is  a  public  domain,  image-processing  and  geographic  informa¬ 
tion  system  originally  developed  by  researchers  in  the  Environmental  Division  of  the  U.S.  Army 
Construction  Engineering  Research  Laboratories  (USACERL)  in  Champaign,  IL.  The  GRASS  sys¬ 
tem  is  used  to  input,  analyze,  and  output  geographic  data  by  users  in  both  military  and  nonmilitary, 
public  and  private  agencies,  based  in  North  America,  Europe,  and  other  parts  of  the  world.  It  is  a 
component  of  the  Integrated  Training  Area  Management  program  implemented  by  USACERL  to 
assist  in  the  management  of  Army  training  lands. 

Objective 

The  objective  of  this  research  was  to  develop  a  map  algebra  language  flexible  enough  to 
accommodate  operations  in  Tomlin’s  first  two  classes  of  map  operations  and  that  would  allow 
developers  and  end-users  to  construct  a  core  of  GIS  analysis  and  image-processing  operations. 

Approach 

A  map  algebra  language  syntax  was  first  defined  and  developed  and  then  tested  and  illus¬ 
trated  with  examples  relevant  to  GISs  and  image-processing  applications.  The  sample  geographic 
data  used  in  the  examples  came  from  the  Spearfish,  SD  database-a  database  distributed  with 
GRASS  software. 

Mode  of  Technology  Transfer 

The  map  algebra  has  been  implemented  in  the  r.mapcalc  module  of  the  Geographic  Re¬ 
sources  Analysis  Support  System.  GRASS  is  transferred  to  the  field  through  the  following  mecha¬ 
nisms:  training  programs,  hands-on  experience,  a  user  documentation  center,  newsletters,  institu¬ 
tional  structures  at  the  Army  and  Interagency  levels,  communication  networks,  and  other  forums. 
GRASS  source  codes  can  be  obtained  via  Internet  computer  software  from  a  file  transfer  protocol 
server  at  ftp  moon.cecer.army.mil.  Current  information  is  available  from  the  GRASS  Information 
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Center  at  USACERL,  which  can  be  contacted  by  phone  at  217-373-7220,  by  fax  at  217-373-7222, 
or  by  e-mail  at  gic@zorro.cecer.army.mil. 


2  LANGUAGE  SYNTAX 

The  r.mapcalc  algebra  employs  mathematical  operators  and  functions  in  expressions  involv¬ 
ing  maps  or  images  that  are  stored  in  raster  grid-cell  format  and  are  two-dimensional  matrices  of 
integer  values. 

The  syntax  for  the  algebra  is  result=expression,  where  expression  is  built  using  maps  and 
images,  mathematical  operators,  functions,  and  temporary  variables.  The  result  map  is  produced 
by  evaluating  the  expression  for  each  cell  in  the  matrix. 

For  example,  the  expression: 


sum  =  mapl  +  map2 

would  produce  a  map  sum  where  each  cell  is  the  sum  of  the  values  of  the  corresponding  cells  in 
mapl  and  map2. 

Maps  and  Images 

Maps  and  images  are  database  files  stored  in  raster  format,  i.e.,  two-dimensional  matrices  of 
integer  values.  The  values  stored  in  maps  may  represent  categorical  (e.g.,  soil  types)  or  continuous 
(e.g.,  elevation)  data,  whereas  the  values  stored  in  images  usually  represent  continuous  data  (e.g., 
satellite  sensor  data).  Although  the  information  represented  by  the  values  within  a  map  or  an 
image  varies,  the  data  format  is  the  same.  Therefore,  this  distinction  between  map  and  image  will 
be  dropped,  and  the  term  map  will  be  used  for  either. 

Map  naming  rules  are  flexible;  however,  a  map  name  must  not  contain  any  special  charac¬ 
ters  (operators,  parentheses,  etc.),  and  it  must  be  distinguishable  from  numbers.  Common  names 
include  elevation,  soils,  vegetation,  roads,  band.l,  etc.,  although  the  format  12oct89  is  also  accept¬ 
able. 

Map  names  may  be  followed  by  a  neighborhood  modifier  which  specifies  a  relative  offset 
from  the  current  cell  being  evaluated.  The  format  is  map[r,c],  where  r  is  the  row  offset  and  c  is  the 
column  offset.  For  example,  map[l,2]  refers  to  the  cell  one  row  below  and  two  columns  to  the  right 
of  the  current  cell,  map[-2,-l]  refers  to  the  cell  two  rows  above  and  one  column  to  the  left  of  the 
current  cell,  and  map[0,l]  refers  to  the  cell  one  column  to  the  right  of  the  current  cell.  This  syntax 
permits  the  development  of  neighborhood-type  filters  within  a  single  map  or  across  multiple  maps. 

At  present,  GRASS  map  files  may  contain  only  integer  values.  However,  GRASS  does 
permit  floating  point  values  to  be  associated  with  each  integer  value  in  a  map.1  Maps  may  be 
prefixed  by  the  @  modifier  which  translates  the  integer  values  in  a  map  to  their  associated  floating 
point  values. 

'These  values  are  stored  in  a  separate  attribute  file  known  as  the  “category  label  file.” 
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Operators 


The  operators  for  the  map  algebra  are: 


operator 

meaning 

type 

precedence 

* 

multiplication 

arithmetic 

4 

/ 

division 

arithmetic 

4 

% 

modulus 

arithmetic 

4 

+ 

addition 

arithmetic 

3 

- 

subtraction 

arithmetic 

3 

== 

equal 

comparison 

2 

i  = 

not  equal 

comparison 

2 

< 

less  than 

comparison 

2 

<  = 

less  than  or  equal 

comparison 

2 

> 

greater  than 

comparison 

2 

>= 

greater  than  or  equal 

comparison 

2 

&& 

and 

logical 

1 

II 

or 

logical 

1 

They  are  applied  from  left  to  right,  with  those  of  higher  precedence  applied  before  those 
of  lower  precedence.  Parentheses  may  be  used  to  control  the  order  of  evaluation.  The  arithmetic 
operations  have  their  usual  meanings,  except  that  division  by  zero  equals  zero.  If  both  operands 
are  integer,  the  result  is  integer,  otherwise  the  result  is  floating  point.  The  comparisons  evaluate 
to  1  if  the  comparison  is  true,  otherwise  they  evaluate  to  0.  The  ||  operation  evaluates  to  1  if 
either  operand  is  nonzero  (true),  and  to  0  otherwise.  The  &&  operation  evaluates  to  1  only  if  both 
operands  are  nonzero  (true),  and  to  0  otherwise. 
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Functions 


Functions  perform  specific  operations  on  a  list  of  expressions  and  result  in  a  single  operand. 
The  r.mapcalc  algebra  includes  the  following  functions: 


abs(x) 

1*1 

exp(x) 

ex 

exp(x,y) 

xv 

log(x) 

In  x 

log(x,y) 

logy* 

sqrt(x) 

atan(x) 

tan'1!  [-90°,90°] 

atan(x.y) 

tan"1!//!  [0°,360°] 

cos(x) 

COS!0 

sin(x) 

sin  !° 

tan(x) 

tanx0 

min(x,y,...) 

minimum  of  {x,y, ...} 

max(x,y,...) 

maximum  of  {x,y, ...} 

float  (x) 

convert  x  to  floating  point 

round(x) 

round  x  to  nearest  integer 

int(x) 

convert  x  to  largest  integer 
less  than  or  equal  to  x 

if(x) 

1,  if  x  !  =  0;  0  otherwise 

if(x,y) 

y,  if  i  !  =  0;  0  otherwise 

if(x,y,z) 

y,  if  x  !  =  0;  z  otherwise 

if(x,p,z,n) 

p,  if  i  >  0;  z,  if  x  ==  0;  n  otherwise 

eval(a,b,...,x) 

x  (but  all  arguments  are  evaluated) 

Temporary  Variables 

Expressions  can  become  complicated  and  values  commuted  in  one  part  of  an  expression 
may  need  to  be  computed  again  in  a  later  part.  Such  a  \  ie  can  be  assigned  a  temporary 
variable  to  be  used  in  place  of  that  value  in  the  later  part  of  the  expression.  For  example. 
result=(map+2)*(map+2)  may  be  expressed,  using  temporary  variable  x,  as  resvlt=(x=map-b2)*x. 
Additional  examples  are  given  on  pages  10  and  13. 


3  GIS  Examples 

The  r.mapcalc  algebra  supports  a  variety  of  GIS  operations,  some  of  which  are  described  in 
the  following  sections;  these,  however,  do  not  exhaust  the  flexibility  and  power  of  the  language.  In 
the  examples,  the  following  maps  are  used: 
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elevation  30-meter  digital  elevation,  range  1066-1840  meters 


rushmore  Camp  Rushmore  (a  fictitious  military  installation) 
0  outside  the  installation 
1  inside  the  installation 


slope  slope  in  degrees 


vegcover  vegetation  cover 


1 

irrigated  agriculture 

2 

rangeland 

3 

coniferous  forest 

4 

deciduous  forest 

5 

mixed  forest 

6 

disturbed 

Map  Recoding 

One  standard  GIS  operation  is  the  recoding  of  values  in  a  map.  For  example,  a  map  of  forest 
cover  versus  nonforest  cover  can  be  created  from  the  vegcover  map  by  recoding  forest  categories  to 
1,  and  nonforest  categories  to  2  as  follows: 


simple.cover  -  if  (vegcover  ==  3  j|  vegcovt  — =  4  ||  vegcover  ==  5,  1)  +  \ 
if  (vegcover  ==  1  ||  vegcover  ==  2  ||  vegcover  ==  6,  2) 

One  adds  the  results  of  both  if  functions  because  each  if  selects  a  distinct  subset  of  the  input  data. 
Note  the  use  of  \  to  indicate  that  the  expression  continues  on  multiple  lines. 

To  make  a  map  of  elevations  between  1200  and  1375  meters  (with  elevations  outside  this 
range  recoded  to  0),  use: 

upland  ss  if  (elevation  >=  1200  &&  elevation  <  =  1375,  elevation) 

Selections 

An  operation  related  to  map  recoding  is  the  selection  or  identification  of  cells  meeting  some 
specified  criteria  in  one  or  more  maps.  Selection  differs  from  recoding  in  two  respects:  (1)  it 
may  involve  more  than  one  map;  and  (2)  the  resulting  map  often  has  only  the  values  0  and  1,  1 
indicating  cells  that  meet  the  criteria  and  0  indicating  those  that  do  not.  The  logical  operators 
(&&  and  jj)  together  with  the  comparison  operators  (==,  !  =,  <,  <  =  ,  >,  and  >=)  provide  this 
selection  capability. 

For  example,  to  create  the  map  upland,  with  elevations  between  1200  and  1375  meters  coded 
as  1,  and  elevations  outside  this  range  coded  as  0,  use: 

upland  =  elevation  >=  1200  &&  elevation  <  =  1375 
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To  create  the  map  forest,  indicating  where  forests  are  found,  use: 


forest  -  vegcover  ==  3  ||  vegcover  ==  4  ||  vegcover  ==  5 


use: 


Then,  to  create  the  combined  map  upland.forest,  indicating  where  upland  forests  aie  found 

upland. forest  =  forest  kk  upland 


Or,  more  directly,  make  the  map  upland.forest  without  first  creating  upland  and  forest  by 

using: 


upland.forest  —  ( elevation  >=  1200  kk  elevation  <=  1375)  kk  \ 

( vegcover  ==  3  ||  vegcover  ==  4  ||  vegcover  ==  5) 

Region  Growing 

A  region-growing  operation  adds  a  one-cell  border  around  an  area  of  a  map  and  can  be 
implemented  with  a  simple  algorithm.  Each  zero  value  in  a  map  is  replaced  by  a  nonzero  value 
from  the  cell  to  the  left,  to  the  right,  above,  or  below. 

For  example,  to  add  a  one-cell  border  to  the  Camp  Rushmore  map,  rushmore ,  one  can 
create  rushmore. grow. 

rushmore.grow  =  if  (rushmore,  rushmore,  \ 

i f(rushmore[0,  -1],  rushmore[ 0,  -1],  \ 

if(rushmore[  0, 1],  rushmore[0 , 1],  \ 

i f(rushmore[—  1, 0],  rus/»more[-l,  0],  \ 
i/(rushmore[l,0],rushmore[l,0]))))) 

The  border  itself,  rushmore. border,  can  now  be  extracted  by  subtracting  the  original  map 
from  the  new  map: 

rushmore. border  =  rushmore.grow  -  rushmore 

Note  that  a  count  of  just  the  border  cells  can  be  used  to  approximate  the  perimeter  of  the 
installation.2 


Slope  and  Aspect 


The  neighborhood  syntax  of  the  r.  mapcalc  algebra  enables  a  user  to  calculate  a  slope  gra¬ 
dient  (the  maximum  rate  of  change  in  altitude)  and  its  aspect  (the  compass  direction  of  the  slope) 
from  elevation  values.  The  basic  formulas  (Dozier  and  Strahler  1983)  are: 


tan(slope) 

tan(aspect) 


\j  (bz/bx)2  +  (bz/iy)2 

6z/6x 

6z/6  y 


2 Cell  connting,  a  regional  operation,  is  not  supported  by  r.mapcalc.  Cell  counts  must  be  determined  using  another 
tool  (e.g.,  the  GRASS  r. stats  command). 
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where  6z/6x  and  6z/6y  are  the  partial  derivatives  in  the  east-west  and  north-south  directions, 
respectively.  Numerical  methods  can  estimate  these  derivatives  (Skidmore  1989).  A  method  given 
by  Horn  (1981)  is: 

\bz/ bx]^^.  =  (zy-liX- 1  +  2ZytX-l  +  Zy+1,1-1  —  Zy-ljc+l  ~  22y,*+l  —  Zy+  l,x+l)  / 8Ax 
[bz/6y\y  x  =  (zy- i,x-i  +  2 Zy-itX  +  2y-iti+i  —  Zy+X'X-i  —  2zy+\tx  ~  ^y+i.i+i )  /8Ay 

where  zy<x  is  the  elevation  value  at  row  y  column  x ,  Ax  is  the  east-west  (i.e.,  column)  grid  spacing, 
and  A y  is  the  north-south  (i.e.,  row)  grid  spacing.  This  can  be  expressed  more  clearly  with  a  matrix 
of  coefficients: 

1  2  1  ' 

0  0  0 
1  -2  -1 
8AX  8  Ay 

The  elevation  map,  with  its  30-by-30-meter  horizontal  grid  spacing,  can  be  processed  with 
r.mapcalc  as  follows: 


slope  =  eval(  x=  (elevation[-l, -1]  +  2  *  elevation[0, -1]  +  e/evatton[l, -1]  \ 

-elevation[- 1, 1]  —  2  *  elevation[ 0, 1]  —  elevation[  1, 1]  \ 

)/(8.0  *  30.0)  ,  \ 

y=  (eIevation[-\,-l]  +  2*elevation[-l,0\  + elevation[-l,l]  \ 
—elevation[  1,  -1]  -2*  elevation[  1,0]  -  elevation[  1, 1]  \ 

)/(8.0  *  30.0)  ,  \ 

atan(sqrt(x  *  x  +  y  *y))  \ 

) 

aspect  =  eval(  x=  (elevation[- 1, -1]  -(-  2  *  elevation[0, -1]  +  elevation[l, -1]  \ 

-elevation[-l,  1]  -  2  *  elevation[0 , 1]  —  elevation[  1, 1]  \ 

)/ (8.0  *  30.0)  ,  \ 

y  =  (elevation[—l,—l]  +  2*elevation[-l,Q]  +  elevation[-l,l]  \ 
—elevation[  1,  -1]  -  2  *  elevation[  1,0]  —  elevation[  1, 1]  \ 

)/(8.0  *  30.0)  ,  \ 

a=  round{atan{x,y))  ,  \ 

if(x\\y,if{a,  a,  360))  \ 

) 


Note  the  use  of  the  temporary  variables  a,  x,  and  y  to  capture  intermediate  results  for  use 
in  a  latter  part  of  the  expression. 

Slope  is  calculated  using  the  single-argument  version  of  atan(),  which  produces  angles  in 
the  range  of  0  to  90  degrees.  For  terrain  with  low  relief,  however,  where  degree  values  may  not 
give  enough  information,  one  could  multiply  the  result  by  10  to  get  lOths  of  a  degree,  or  omit  the 
atan( )  function  to  get  tan(slope )  instead. 

Aspect  is  calculated  using  the  two-argument  version  of  atan(),  which  produces  angles  in  the 
range  of  0  to  360  degrees.  However,  because  aspect  is  undefined  if  both  x  and  y  are  zero  (i.e.,  flat 
terrain),  additional  care  is  required.  The  code: 

a  =  round(atan(x ,  y)) 

*f(x\\y,  if  (a,  a,  360)) 
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produces  values  from  1  to  360  for  cells  that  have  nonzero  slope,  and  a  zero  value  for  flat  terrain. 
Frank  (1988)  used  an  alternative  technique  to  compute  hzjbx  and  6z/6y: 


6z/6x  = 


0  0  0 
0  0  0 
1  8  0 
0  0  0 
0  0  0 


12AX 


6z/6y  = 


0  0  10  0 

0  0  8  0  0 

0  0  0  0  0 

0  0-100 
00-800 
12A  Y 


This  is  expressed  in  r.mapcalc  as: 

slope  =  eval(  x  =  ( 


x  =  (elevation[ 0,  —2]  +  8  *  elevation[ 0,  —1] 
—elevation[0, 2]  —  8  *  elevation[0, 1] 

)/ (12.0  *  30.0)  , 

y  =  (elevation[— 2,0]  +  8*  e/euatton[-l,0] 
-elevation[2, 0]  —  8  *  elevation[l,Q\ 
)/( 12.0*30.0), 
atan(sqrt(x  *  x  +  y  *y)) 


aspect  =  eval(  x  = 


x  =  ( elevation[0 ,  —2]  +  8  *  elevation[ 0,-1] 

-e/evatton(0,2]  —  8  *  elevation[Q,  1] 
)/(12.0*  30.0)  , 

y  =  (elevation[— 2, 0]  +  8  *  elevation[— 1, 0] 
-elevation[ 2, 0]  -  8  *  elevation[l,  0] 
)/(12.0*  30.0)  , 
a=  round(atan(x,  y))  , 
if(x\\y,  if  (a,  a,  360)) 


Hydrologic  Simulation 

r.mapcalc  can  construct  a  simple  hydrologic  model  that  iteratively  “flows”  water  across  a 
landscape.  A  constant  amount  of  water  is  first  deposited  in  each  cell  in  the  landscape.  Then,  at 
each  time-step,  a  portion  of  the  water  in  each  cell  is  drained  into  its  eight  surrounding  neighbors. 
The  basic  logic  is: 


for  each  cell 

for  each  neighbor 

let  height  difference  =  (cell  elevation  +  water  height) 
minus  (neighbor  elevation  +  water  height) 
if  the  height  difference  is  positive 
then 

if  the  cell  elevation  is  greater  than 

the  neighbor’s  elevation  +  water  height 

then 
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drain  out  a  portion  of  the  water 
otherwise 

drain  out  a  portion  of  the  height  difference 

otherwise 

if  the  neighbor's  elevation  is  greater  than 
the  cell's  elevation  +  water  height 

then 

drain  in  a  portion  of  the  neighbor’s  water 
otherwise 

drain  in  a  portion  of  the  height  difference 

end 

end 


The  outer  loop  is  for  each  cell  in  the  landscape  and  is  done  automatically  by  r.mapcalc. 
The  inner  loop  applies  to  the  eight  neighboring  cells  and  must  be  explicitly  coded  by  the  user. 

Before  running  this  model,  the  elevation  map  must  be  converted  to  the  same  units  as  the 
water  map  and  must  be  filtered  to  smooth  the  elevation.  In  this  example,  let  us  assume  that  the 
water  will  be  in  units  of  0.1  inches.  To  convert  elevation  from  meters  to  0.1-inch  units: 

elev  =  elevation  *  393.7 


The  map  elev  can  then  be  smoothed  with  the  following  filter: 


1  2  1 
2  8  2 
1  2  1 

20 


as  follows: 


elev  =  (  elev[-l, -1]  +  2  *  e/en[-l,0]  -I-  e/er[-l,  1]  \ 

+  2  *  elev[0,-l]  + 8*  elev +  2  *  elev[0,l)  \ 
+  e/ev[l, -1]  +  2  *  e/ev[l,0]  +  elev[l,  1]  \ 

)  /20 


(This  filter  produces  a  map  with  smoother  contours.  It  also  tends  to  create  dams  at  terrain  choke 
points,  which  will  form  little  lakes  during  the  simulation.) 
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The  model  itself  is  coded  as  follows: 


water  =  water  +  eval(x  =  elev  +  water ,  \ 


if 

(x>(y  = 

e/etr[ — 1, 0]  +  water[ 

-i,o]), 

\ 

-.15  *if 

(elev  >  y, 

water , 

x  -  y), 

\ 

.15*  if 

(e/eu[-l,0]  >  x, 

water[— 1,0], 

y-x))+ 

\ 

if 

(*>(»  = 

e/ev[l,0]  -1-  water[l,0]), 

\ 

-.15*t/ 

(elev  >  y, 

water , 

x-y). 

\ 

.15**/ 

(e/ev[l,0]  >  x, 

water  [1,0], 

y-x))+ 

\ 

if 

(x>(y  = 

e/ev[0,  -1]  +  water[0,  -1]), 

\ 

-.15**/ 

(elev  >  y, 

water , 

x-y), 

\ 

.15*  if 

(eiet>[0,  -1]  >  x, 

water[0,  -1], 

y-x))+ 

\ 

if 

(x>(y  = 

e/et>[0, 1]  -1-  water[Q, 1]), 

\ 

-.15  *  if 

(elev  >  y, 

water, 

x-y), 

\ 

.15  *  if 

(elev[ 0, 1]  >  x, 

water[  0, 1], 

y  -  x))+ 

\ 

if 

(a?  >  (y  = 

ele  ®[— 1, 1]  +  water[- 

-1,1]), 

\ 

— .10  *  */ 

(elev  >  y , 

water, 

x-y), 

\ 

.10**/ 

(e/ev[-l,  1]  >  x, 

water[—  1, 1], 

y-x))+ 

\ 

if 

(x  >  (y  = 

ele t>[l,  1]  -|-  itfafer[l,  1]), 

\ 

-.10**7 

(elev  >  y, 

water, 

x-y), 

\ 

.10**7 

(e/et>[l,  1]  >  x, 

water[l,  1], 

y-x))+ 

\ 

if 

(*>(»  = 

elev[l,  -1]  +  water[l ,  -1]), 

\ 

-.10**/ 

(elev  >  y, 

water, 

x-y), 

\ 

.10**/ 

(e/et)[l,  -1]  >  x, 

water  [1,-1], 

y  -  x))+ 

\ 

if 

(x>(y  = 

e/et>[-l,  -1]  +  ti?ater[-l,  -1]), 

\ 

—.10  *  if 

(elev  >  y, 

water. 

x-y), 

\ 

.10  *  if 

(e/e«[-l,-l]  >  x, 

water[-l,  -1], 

y-x))) 

The  resulting  map  water  represents  the  water  depth  after  a  single  iteration  (time-step). 
Note  that  the  eight  sections  in  this  model  (beginning  with  t f(x  >  . . .))  are  similar,  each  handling 
a  neighboring  cell.  The  first  four  drain  15  percent  of  the  water  to  or  from  the  cells  above,  below, 
and  to  either  side;  the  last  four  drain  10  percent  of  the  water  to  or  from  the  cells  at  the  diagonal 
positions. 

This  model  has  been  applied  to  a  section  of  the  Spearfish  database  with  realistic-looking 
results.  (The  data  set  entitled  spearfish  is  a  sample  database  included  with  the  GRASS  release 
and  covers  the  Spearfish,  South  Dakota  vicinity.)  Water  drains  away  from  the  uplands  forming 
temporary  streams  and  ponds.  For  those  wishing  to  repeat  the  experiment  with  GRASS,  enter  the 
model  code  in  a  file  called  water.mapcalc,  then  create  a  controlling  shell  script: 

#!/bin/sh 

r.mapcalc  water=120  #  start  with  12  inches  of  waterin  each  cell 

d.rast  water  #  display  the  water  map 

i=l 

while  [  $i  !=  100] 
do 


n=l 

while  [  $n  !=  10  ] 
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do 


r.mapcalc  <  water. mapcaJc  #  run  the  simulation 
d.rast  water  #  display  the  water  map 

n =‘expr  $d  +  1‘ 

done 

g.copy  ras t = water, water. $i  #  save  a  snapshot 

i=‘expr  $i  +  V 

done 

This  script  will  run  1000  iterations.  The  depth  starts  at  120  (12  inches,  or  30.48  cm,  of 
water)  across  the  entire  map.  After  every  tenth  iteration,  the  resultant  map  is  copied  to  a  map 
with  a  unique  name3  for  rapid  display  later.  One  can  run  the  shell  script  and  watch  the  water 
“flow”  off  the  slopes  of  South  Dakota’s  Black  Hills. 

This  model  could  be  improved  by  considering  additional  factors  such  as:  (1)  evaporation  and 
transpiration — in  each  time-step  (iteration),  some  water  will  evaporate,  which  could  be  modeled 
either  as  a  constant  value  or  as  a  simple  function  of  groundcover  type;  (2)  flow  impedance — the 
amount  of  water  that  can  leave  a  cell  in  a  given  time  step  could  be  modified  using  groundcover 
information;  and  (3)  ground  saturation— the  rate  at  which  soils  absorb  water  during  rainfall  can  be 
estimated  based  on  soil  type  and  length  of  exposure  to  water. 

4  IMAGERY  EXAMPLES 

The  r.mapcalc  algebra  also  supports  many  image-processing  operations.  Examples  of  these 
follow;  again,  they  are  not  exhaustive,  but  are  intended  to  illustrate  the  power  of  the  algebra. 

Spectral  Ratios 

Ratio  transformation  of  spectral  data  is  one  technique  used  in  the  analysis  of  remotely  sensed 
data  (Lillesand  and  Kiefer  1987).  Between-band  ratios,  which  eliminate  multiplicative  effects,  are 
easily  created  with  r.mapcalc: 

ratio  =  band.l/band.2 

Ratios  of  between-band  differences,  which  eliminate  additive  effects,  are  also  straightforward: 

ratio  =  (band.  1  —  band. 2) / (band. 3  -  band. 2) 

A  normalized  vegetation  index  ( nvi )  for  AVHRR  data  is  computed  as: 

nvi  =  (avhrr.2  —  avhrr.l)/(avhrr.2  +  avhrr.l) 

These  results  could  be  produced  with  a  very  simple  algebra  that  uses  only  arithmetic  operators. 
The  transformed  vegetative  index  ( tvi)  for  the  Landsat  Thematic  Mapper  illustrates  the  use  of 
functions  in  r.mapcalc: 


tvi  =  100  *  sqrt((tmA  -  tm.3)/(tm.4  +  tm. 3)  +  0.5) 

3These  maps  will  be  named  water.  1,  water. 2,  . . .  water.  100. 
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Note  that  0.5  is  added  to  help  ensure  that  sqrtQ  is  not  applied  to  a  negative  value.  This  is 
guaranteed  by  using  the  max()  function: 

tvi  =  100  *  sqrt(max( 0,  (fm.4  -  tm.3)/(tmA  +  tm. 3)  +  0.5)) 

If  the  results  are  to  be  in  the  range  0  to  255,  this  ratio  will  produce  such  results: 

ratio  =  255.0/90.0*  atan(float(tm.  1)/ float{tm. 2)) 

The  float()  function  ensures  that  the  division  is  floating  point,  not  integer,  division.  Because  the 
output  map  can  hold  only  integer  values,  r.mapcalc  will  round  the  results  to  integers  before  they 
are  stored  in  the  map. 

Spatial  Filters 

A  common  image-processing  operation  is  local-neighborhood  filtering.  A  small  window  is 
moved  over  an  image,  and  the  center  pixel  is  replaced  by  a  combination  of  all  the  pixels  in  the 
window. 

Low-frequency  filters  (also  called  low-pass  filters)  deemphasize  high  spatial  frequencies  and 
involve  computing  a  weighted  average  of  all  pixels  in  the  window.  A  simple  average,  which  smoothes 
the  image,  is  based  on  the  filter: 


9 

This  is  expressed  with  r.mapcalc  as: 

ave  3x3  =  (  image[— 1, -1]  +  image[— 1,0]  +  image[-l,  1]  \ 

+  image[ 0,  -1]  +  image[ 0, 0]  +  image[ 0, 1]  \ 

+  tmapejl,— 1]+  »ma(/e[l,0]  +  imu^e[l,l]  \ 

)  /  9 

Use  high-frequency  filters  (high-pass  filters)  to  emphasize  high  spatial  frequencies  and  for 
edge  enhancement.  The  technique  is  the  same  as  for  low-frequency  filters,  except  that  some  of  the 
weights  are  negative.  One  such  filter  is: 

-1  -1  -1 
-1  9  -1 

-1  -1  -1 

This  is  expressed  with  r.mapcalc  as: 

high3x3  =  (  -image[- 1,—  1]  —  image[— 1,0]  -  image[- 1,1]  \ 

-image[ 0,  - 1]  +  9  *  image[ 0, 0]  -  image[ 0, 1]  \ 

-image[  1,  -1]  -  ima</e[l,0]  -  image[  1, 1]  \ 

) 
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The  Sobel  filter  is  a  nonlinear  edge-enhancement  filter: 


-1  0  1 

-2  0  2 

-1  0  1 


1  2  1 

0  0  0 

-1  -2  -1 


It  is  similar  to  the  method  to  determine  slope  values  from  elevation  values  and  is  expressed  in 
r.mapcalc  as: 

sobelSxS  =  eval(  x  =  image[- 1, 1]  +  2  *  image[ 0, 1]  -I-  image[  1, 1]  \ 

-ima<7e[-l,  —  1]  -  2  *  image[0,  -1]  -  image\\,  —1]  ,  \ 
y  =  imaffe[-l, -1]  +  2*  *masre[-l,0]  +  tma0e[-l,l]  \ 

— image[l,  -1]  -  2  *  »ma5e[l,0]  —  tmage[\,  1]  ,  \ 

sqrt(x  *x  +  y*y)  \ 


Another  type  of  filter  is  an  adaptive  one,  which  replaces  the  center  pixel  only  if  some 
criteria  are  met.  Eliason  and  McEwen  (1990)  describe  two  whereby  the  center  pixel  is  replaced  by 
the  average  value  in  the  neighborhood  if  the  difference  between  them  exceeds  a  threshold.  One 
filter  sets  the  threshold  to  a  multiple  of  the  statistical  variance  in  the  neighborhood.  The  other 
sets  the  threshold  to  a  multiple  of  the  statistical  variance  in  the  entire  image.  Both  filters  can  be 
coded  with  r.mapcalc ;  the  variance  within  the  entire  image  is  a  regional  calculation  and  must  first 
be  computed  by  another  command4  and  then  inserted  into  the  r.mapcalc  code.  The  first  filter  is 
illustrated. 

Let  P  be  the  value  of  the  center  pixel,  s  the  sum  of  all  values  in  the  neighborhood,  n 
the  number  of  pixels  in  the  neighborhood,  and  ss  the  sum  of  the  squares  of  the  values  in  the 
neighborhood.  The  average  is  given  by: 


and  the  variance  is  given  by: 


The  center  pixel  is  replaced  if: 


s 

P  =  ~ 


o  SS 

"  =n~' 


(P  -  p?  >  Ccr2 


where  C  is  usually  between  1.0  and  4.0.  The  r.mapcalc  code  for  a  3x3  window  with  C  set  to  2.25 


4The  GRASS  r.stat$  command  could  be  used  to  compute  the  variance  across  the  entire  image. 
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filter  =  eval(  s  =  image[-\, -1]  +  ima0e[-l,O]  +  tmagre[-l,  1]+ 

image[0,  -1]  +  tma^e[0, 0]  +  image[ 0, 1]+ 
image[\,  -1]  +  tmage[l,0]  +  image[\,  1]  , 
ss  -  image[—  1, -1]  *  image[—  1,  — 1]+ 
image[- 1,0]  *  image[- 1,0]+ 
image[- 1, 1]  *  tmage[- 1, 1]+ 
tmajfe[0,  -1]  *  image[ 0,  -1]+ 
tma<7e[0,0]*  image[ 0,0]+ 
image[ 0, 1]  *  image[ 0,  l]+ 
image[  1,  —1]  *  image[\,  — 1]+ 
image[l,0]  *  irna<7e[l,0]+ 
image[l,  1]  *  image[l,  1]  , 
ave  =  s/9.0  , 
var  =  ss/ 9.0  -  ave  *  ave  , 
x  =  image  —  ave  , 
if(x  *  x  >  2.25  *  var ,  ave ,  image) 

) 


\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 


Radiometric  Calibration 

Radiometric  calibration  converts  digital  values  recorded  by  a  remote  sensing  system  to 
radiance  values  at  the  top  of  the  atmosphere,  possibly  followed  by  a  conversion  to  surface  reflectance 
values.  Radiance  values  can  be  computed  using  the  following  formula  (Chavez  1989): 

radi(x,y)  =  (dn,(x,  y)  -  offseti)/gairii 

where  rad  is  the  radiance  for  band  i  at  pixel  x,y,  dn  is  the  value  recorded  by  the  sensors;  and  offset 
and  gain  are  the  sensor  offset  and  gain  for  band  i.  This  formula  is  easily  handled  by  r.mapcalc. 
For  example,  using  the  values  for  the  offset  and  gain  as  reported  by  Chavez  for  October  3,  1988 
Thematic  Mapper  data,  radiance  images  are  constructed  as  follows: 

retd.l  =  (tm.  1  -  2.4899)/ 16.5993 
rad. 2  =  (fm.2  -  2.3871)/8.5104 
rad. 3  =  (tm. 3  -  1.4815)/12.4074 
radA  =  (tm.4  -  1.8418)/12.2790 
rad.b  =  (tm. 5  -  3.4240)/92.5292 
rad. 7  =  (tm.7  -  2.6323)/175.4878 

The  radiance  formula  for  the  SPOT  sensor  appears  in  a  slightly  modified  form  (USA  SPOT  Image 
Corporation  1989): 

radi(x,y )  =  dm(x,y)/Ai 

where  Ai  are  the  absolute  calibration  coefficients.  The  May  27,  1989  SPOT  multispectral  image  in 
the  Spearfish  database  can  be  converted  to  radiance,  using  r.mapcalc ,  as  follows: 

rad.l  =  5potm5.1/1.05586 
rad.2  =  spot,  ms.2/1. 12140 
rad.3  =  spot.ms. 3/0.97244 
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Surface  reflectance  calculations  are  more  involved  because  the  formulas  must  incorporate 
atmospheric  effects.  There  are  various  methods  for  modeling  these  effects  (Richards  1986,  Price 
1987,  Chavez  1989);  a  formula  given  by  Chavez  is: 

R  (x  )  -  Kdist2[radi{x,y)  -  feoze,) 

1  ’  E{slope(x ,  y)  ■  sun  ■  mhaze , 


where  R  is  the  surface  reflectance  image,  rad  is  the  radiance  image,  slope  is  a  slope  map,  dist  is  the 
Earth-Sun  distance,  haze  and  mhaze  are  the  additive  and  multiplicative  atmospheric  haze  factors, 
E  is  exoat  mo?  nheric  spectral  irradiance,  and  sun  is  the  sun  elevation  at  the  time  the  image  was 
captured.  C<  ing  the  nonmap  factors  into  single  constants,  the  formula  can  be  written  in  the 


form: 


Ri{x,y) 


Aj[radj(x,  y)  -  ifrj 
slope{x,y) 


which,  after  substituting  the  appropriate  values  for  the  constants  A  and  5,  is  expressed  with 
r.mapcalc  as: 


TCI  =  A\  *(rad.l- B\)/ slope 
R.2  =  A2*  (rad.  2  -  B^)  I  slope 
R. 3  =  A3  *  (rad.3  -  B3)/ slope 


Principal  Components 

Principal  components  analysis  involves  transforming  a  set  of  correlated  variables  into  a  new, 
uncorrelated  set.  The  variables  are  the  band  files  from  a  multispectral  sensor.  The  transformation 
is  a  linear  combination  of  the  band  data: 

N 

component  =  *  band]  i  —  1, . . . ,  N 

i=i 

where  N  is  the  number  of  band  files  and  tn,j  are  the  eigenvectors  for  the  between-band  covariance 
matrix.  Although  neither  the  covariance  computation  nor  the  determination  of  the  eigenvectors  can 
be  accomplished  using  the  map  algebra,  the  components  can  be  computed.  For  example,  suppose 
the  three  multispectral  bands  from  a  SPOT  image  had  the  following  covariance  matrix: 

'  462.88  480.41  281.76  ' 

480.41  513.02  278.92 
281.76  278.91  336.33 

A  set  of  eigenvectors  for  this  matrix  (in  decreasing  order  of  spectral  variance)  is: 

vectorj  21.24  22.15  14.77 

vector  2.91  4.46  -10.87 

vector  1.82  -1.62  -0.18 

To  compute  the  second  component  with  r.mapcalc: 

pc. 2  =  2.91  *  spot.ms.l  +  4.46  *  spot.ms.2  -  10.87  *  spot.ms. 3 
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Merging  Panchromatic  Data  With  Multispectral  Data 

The  French  SPOT  satellite  produces  a  10-meter  panchromatic  image  plus  three  spectral 
bands  with  20-meter  resolution.  Landsat  Thematic  Mapper  produces  six  spectral  bands  with  30- 
meter  resolution  (as  well  as  a  120-meter  resolution  thermal  band).  A  sharpened  color  image  can 
be  produced  by  merging  the  high- resolution  panchromatic  image  with  the  lower  resolution  spectral 
bands.  With  the  assumption  that  the  spectral  data  have  been  registered  and  resampled  to  the 
10-meter  panchromatic  data,  a  number  of  techniques  have  been  used  to  perform  the  merger,  most 
of  which  can  be  done  with  the  map  algebra. 

Welch  and  Ehlers  (1987)  described  two  direct  methods: 

M[  =  di(wiMi  +  W2P)  +  bi 
Ml  =  at{MtP)$  +bt 

is  spectral  band  i,  P  is  the  panchromatic  band,  w\  and  W2  are  weights,  and  a,  and  b,  are 
make  the  results  fall  within  the  range  0  to  255.  Both  these  methods  can  be  implemented 
r.mapcalc  algebra. 

simple  average,  based  on  the  first  method,  is  coded  as  follows: 

merge.  1  =  (spot.ms.  1  +  spot.pan)/ 2 
merge.  2  =  (spot.  ms.  2  +  spot.pan)/ 2 
merge.Z  =  (spot.ms. 3  -I-  spot.pan)/ 2 

Or  a  combination  of  both  methods  can  be  used: 

merge.  1  =  sqr t(spot. pan  *  spot.ms.  1) 
merge.2  =  sqrt(spot.pan  *  spot.ms.2) 
merge. 3  =  (spot.pan  +  3*  spot. ms. 3)/ 4 

A  third  method  transforms  three  selected  bands  from  red,  green,  and  blue  color  space  to 
intensity,  hue,  and  saturation  color  space.  The  intensity  is  replaced  by  the  panchromatic  image, 
and  the  new  intensity,  hue,  and  saturation  combination  is  transformed  back  into  red,  green,  and 
blue. 

Intensity,  Hue,  Saturation 

Conversion  from  red,  green,  and  blue  colors  to  intensity,  hue,  and  saturation  involves  for¬ 
mulas  that  can  be  represented  using  the  r.mapcalc  algebra.  The  following  is  based  on  the  hue, 
saturation,  and  value  (HSV)  algorithm  found  in  Foley  and  Van  Dam  (1984).  Assuming  that  R,  G, 
and  B  represent  the  red,  green,  and  blue  components  of  an  image,  the  value  (V'),  saturation  (5), 


where  Mi 
chosen  to 
using  the 

A 
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and  hue  (H)  can  be  formulated  as  follows: 


V  =  max(R,G,B) 

5  =  255.0  *(V  -min(R,G,B))/V 


H  =  if(S,eval(  \ 

m  =  float(min(R,G,  B))  ,  \ 

r  =  (V  -  R)/(V  -  m)  ,  \ 

g  =  (V  -  G)/(V  -  m)  ,  \ 

b  =  (V-B)l(V-m),  \ 

h  =  if(R  ==  V,b-  g)+  \ 


if(G==V,2+r-b)+  \ 
if(B==  V,4  +  5-r),  \ 
if{h  <=  Q,h  +  6,h)  +  60)) 

Although  saturation  is  normally  defined  in  the  range  from  0  to  1,  it  is  multiplied  here  by  255  to 
retain  more  information  after  conversion  to  integer.  Hue  will  be  in  the  range  of  0  to  360,  with  0 
representing  undefined  hues  (i.e.,  white,  black,  or  gray). 

The  inverse  transformation  is: 


R  —  eval(  \ 

s  =  5/255.0,  \ 

h  =  tf(H  >=  360,  H  -  360,  if  )/60.0  ,  \ 

i  =  int(h)  ,  \ 

f  =  h-i  ,  \ 

P  =  y*(l-s),  \ 

q  =  V*(l  -s*f),  \ 

t  =  V  *{l- s*(l-  f))  ,  \ 

*/( *  ==  0,  V)  +  if(i  ==  1,9)  +  »/(»'  ==  2,p)+  \ 
*'/(*'  ==  3 ,p)  +  if(i  =-  4,  t)  +  if(i  ==  5,  V)) 

G  =  eval(  \ 

s  =  5/255.0  ,  \ 

h  =  if(H  >=  360,  H  -  360,  tf)/60.0  ,  \ 

t  =  int(h)  ,  \ 

/  =  h-i,  \ 

P  =  V*(l-s)  ,  \ 

q  =  v*(i-s*f),  \ 

t  =  v  *(1  —  s  +  (l  —  /)) ,  \ 

if(i  ==  0,0  +  */(*  ==  1,V)  +  */(«  ==  2,V)+  \ 
if(i  ==  3  ,q)  +  t/(t  ==  4,p)  +  if(i  ==  5,p)) 
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B  =  eval(  \ 

s  =  5/255.0,  \ 

h  =  if(H  >=  360, /?  -  360, //)/60.0  ,  \ 

i  =  int(h)  ,  \ 

/  =  h-i ,  \ 

p  =  V  *  (1  -  s)  ,  \ 

q  =  V*(l-a*f),  \ 

t  =  K*(l-s*(l-/)),  \ 


if(i  ==  o  ,p)  + 1/(»  ==  i,p)+  */(*  ==  2>0+  \ 
if(i  ==  3,  V)  +  i/(t  ==  4,  V)  +  */(*  ==  5,q)) 


5  CONCLUSION 

Providing  a  map  algebra  for  manipulation  of  raster  map  data  results  in  a  flexible  geographic 
information  system.  Given  appropriate  raster  data  layers,  the  type  and  number  of  potential  data 
manipulations  are  virtually  limitless.  For  many  applications,  users  are  freed  from  software  limits 
and  are  bound  only  by  their  ability  to  employ  the  algebra,  r.mapcalc  supports  three  levels  of  usage: 
(1)  as  a  resource  for  users  who  need  to  perform  specific  algebraic  functions  that  are  not  provided 
by  other  GRASS  programs;  (2)  as  a  foundational  tool  for  advanced  GIS  users  to  develop  a  limitless 
set  of  image  and  map  analysis  functions;  and  (3)  as  a  resource  for  programmers  and  developers  to 
use  in  building  new  functions. 

r.mapcalc  is  perhaps  most  important  at  the  intermediate  level,  as  a  tool  for  advanced  GIS 
users.  A  common  implementation  model  for  a  GIS  in  larger  organizations  is  one  or  two  sophisticated 
users  supporting  numerous  casual  users.  In  this  context,  r.mapcalc  becomes  a  language  that  opens 
an  array  of  possibilities  for  spatial  analysis  and  provides  a  means  to  develop  macros  to  support 
routine  operations  used  by  less  sophisticated  users. 
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